#include "LocalGeneralizer.h"
#include "KNode.h"

using namespace NEWMAT; 



virtual double LocalAverager::addData( KData & pt)
{
	if (pt.second > max)
		max = pt.second; 

	if (pt.second < min)
		min = pt.second; 

	mean = ((mean*size) + pt.second )/ (size+1); 
	size++; 
	return computeError(); 

}

double LocalAverager::computeError()
{
	double err = fabs(max - mean); 
	if (fabs(min - mean) > err)
		err = fabs(min-mean); 

	return err; 

}


virtual double LocalAverager::addData(list<KData> &pts)
{
	double sum = mean*size; 
	size += pts.size(); 

	for(list<KData>::const_iterator it = pts.begin(); it != pts.end(); it++)
	{
		if ((*it).second > max)
			max = (*it).second; 

		if ((*it).second < min)
			min = (*it).second; 

		sum += (*it).second ; 
	}
	mean = sum/pts.size(); 

	return computeError(); 
}




void LinearSolver::updateWeights()
{
   Matrix U, V; DiagonalMatrix D;

   SVD(X,D,U,V);                              // X = U * D * V.t()

	W= 0.0; 
	
//	printf(" eigen values\n"); 
//	cout << setw(4) << setprecision(6) << D; 

	for(int i=1; i<= dim+1 ; i++)
		if ( D.element(i-1) > 0.01) 
			W = W + (dotproduct(U.column(i), Y)/ D.element(i-1))*V.column(i); 

	
}

LinearSolver::~LinearSolver()
{
#ifdef DEBUG_PRINT
	printf("deleting linearsolver with size %d\n", size); 
#endif

//	cout << setw(4) << setprecision(6) << (X | Y )  ; 
//	cout << setw(4) << setprecision(6) << W; 

}


double LinearSolver::addData(KData &pt)
{
	size++; 

	if (size == 1)
	{
		X.resize(size, dim+1); 
		Y.resize(size); 
	}else 
	{
		X.resize_keep(size, dim+1); 
		Y.resize_keep(size); 
	}

	X.row(size) << pt.first ; 
	X.element(size-1,dim) = 1.0; 
	Y.element(size-1) = pt.second ; 

	updateWeights(); 
	return getLSError(); 
}

double LinearSolver::addData(list<KData>& pts)
{
	size += (int) pts.size(); 

	if (size == (int)pts.size())
	{
		X = Matrix(size, dim+1);	//the last column is always 1
		Y = ColumnVector( size); 
	}else
	{
		X.resize_keep(size, dim+1);
		Y.resize_keep(size); 
	}

	int i=size - (int) pts.size()+1; 
	for(list<KData>::const_iterator  it = pts.begin(); it!= pts.end(); it++)
	{
		X.row(i) << (*it).first ;
		X.element(i-1, dim) = 1.0;		//this is always 1
		Y.element(i-1) = (*it).second ; 
		i++; 
	}

	updateWeights(); 

	Matrix res = (X*W); 

#ifdef DEBUG_PRINT
	cout << setw(dim) << setprecision(3) << (W ) ; 
	cout << setw(4) << setprecision(6) << (X  |  Y | res) ; 
#endif

	return getLSError(); 
}


double LinearSolver::getLSError()
{
	return (X*W - Y).norm_infinity(); 
}

